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Abstract. The Density Matrix Renormalization Group (DMRG) was introduced by 
Steven White in 1992 as a method for accurately describing the properties of one- 
dimensional quantum lattices. The method, as originally introduced, was based on 
the iterative inclusion of sites on a real-space lattice. Based on its enormous success 
in that domain, it was subsequently proposed that the DMRG could be modified for 
use on finite Fermi systems, through the replacement of real-space lattice sites by 
an appropriately ordered set of single-particle levels. Since then, there has been an 
enormous amount of work on the subject, ranging from efforts to clarify the optimal 
means of implementing the algorithm to extensive applications in a variety of fields. In 
this article, we review these recent developments. Following a description of the real- 
space DMRG method, we discuss the key steps that were undertaken to modify it for 
use on finite Fermi systems and then describe its applications to Quantum Chemistry, 
ultrasmall superconducting grains, finite nuclei and two-dimensional electron systems. 
We also describe a recent development which permits symmetries to be taken into 
account consistently throughout the DMRG algorithm. We close with an outlook for 
future applications of the method. 
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1. Introduction 

The properties of non-relativistic many-body quantum systems, as arise in many 
branches of physics, are described by the Schrodinger equation. Unfortunately, the 
Schrodinger equation can be solved exactly in but a few exceptional cases, either 
involving one- dimensional or symmetry-restricted Hamiltonians or for systems with 
very few particles. For almost all problems, reliable approximate solutions must be 
sought. Many methods of approximately solving the Schrodinger equation for many- 
body quantum systems have been developed, often tailored to the specific system or 
problem under investigation. Especially challenging has been to come up with methods 
that can treat large numbers of particles with great accuracy, particularly in the vicinity 
of a phase transition. A major breakthrough was the development of the numerical 
renormalization group (RG) method by Wilson pQ in the mid-1980s, which provided for 
the first time an efficient strategy for solving the Kondo problem The enormous 
power of the RG method immediately suggested its application to quantum lattices. 
It soon became clear, however, that the Wilson RG method ran into difficulties when 
dealing with such systems. 

The Density Matrix Renormalization Group (DMRG) was introduced by White 
in 1992 [3 EJ E| in an effort to overcome the limitations of Wilson's RG in describing 
one-dimensional (ID) quantum lattice models. The new method was soon shown to be 
extremely powerful, producing results for the ground state energy of the S=l Heisenberg 
chain that were accurate to twelve significant figures pE], well beyond the precision of 
large-scale diagonalization methods combined with finite size corrections or of Monte 
Carlo techniques. 

The DMRG is an approximate variational procedure that is rooted in Wilson's onion 
picture jSj. In the context of quantum lattices, the idea is to start with a set of lattice 
sites and then to iteratively add to it subsequent sites until all have been treated. At each 
iteration, the dimension of the enlarged space increases as the product of the dimension 
of the initial subspace and that of the added site. The renormalization group procedure 
consists of reducing the enlarged space to the same dimension as the initial subspace 
and then transforming all operators to this new truncated basis (renormalization). 

One of the key features that distinguishes the DMRG from Wilson's RG is the 
criterion by which the truncation is implemented. While the Wilson RG retains the 
lowest Hamiltonian eigenstates of the enlarged space, the DMRG uses a very different 
strategy. Here the idea is to construct the reduced density matrix for the enlarged space 
in the presence of a medium that approximates the rest of the Hilbert space, then to 
diagonalize this density matrix, and finally to maintain only those states with the largest 
density matrix eigenvalues. This method of truncation is guaranteed to be optimal in 
the sense that it maximizes the overlap of the approximate (truncated) wave function 
with the wave function prior to truncation. 

Following the ideas just described, one can treat the entire set of "onion" layers. 
Depending on the manner in which correlations between layers fall off and in the 
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choice of the order in which layers are treated, this can sometimes lead to an accurate 
representation of the ground state of the system. Usually, it does not, however, since the 
early layers know nothing of the physics of those treated subsequently. This suggests the 
use of a "sweeping" algorithm, whereby once all layers have been sampled, we simply 
reverse direction and update them based on the information of the previous "sweep". 
Such a sweeping algorithm can be iteratively implemented until acceptable convergence 
in the results has been achieved. 

The real-space DMRG (rDMRG), as described above, has been applied extensively 
and with enormous success to many ID systems 0|7j, such as spin chains and ladders 
[HI Ej, t-J [TO] and Hubbard models [TT], models with disorder [12], models with 
impurities ^3], polymers [H], and the Kondo model [T^]. Extension of the method 
to 2D lattices ^H] requires establishing a path that covers all of the sites of the lattice. 
Since there are many possible paths through a 2D lattice, it is critical that the path 
chosen conform to the dominant correlations of the problem. Electrons in small square 
or rectangular lattices with t-J or Hubbard Hamiltonians have been treated using 
this DMRG strategy, though with much less precision than for ID systems. Problems 
in 3D are beyond the scope of the DMRG in real space. 

A pedagogical introduction to the real space DMRG method together with a 
summary of some of its most important applications can be found in the proceedings 
of the workshop "Density-Matrix Renormalization" ^H] held in Dresden in September 
1998. 

Since then, there have been many important developments and advances of the 
DMRG, which has opened up the possibility of its application to finite Fermi systems. 
These new developments originated with the introduction of a transformation of the 
DMRG from real space to momentum space (kDMRG). The kDMRG method was first 
discussed by Xiang JH] and applied to the Hubbard model in 2D. The results were not 
very encouraging, however, and the method was therefore not pursued at that time for 
other 2D quantum lattices. 

This method was recovered in a slightly modified form in 1999 in an effort to build 
a DMRG methodology appropriate to problems in quantum chemistry [20, 21 J . In this 
approach, the momentum basis was replaced by a basis of single-particle orbitals, so 
as to better accommodate the dynamics of the problem. The first applications were 
reported in refs. [211] 121 j . Several other applications to problems in quantum chemistry 
were subsequently reported [221 EH] 

A somewhat different extension of the DMRG methodology was proposed in refs. 
|24| 125] , in the context of a study of ultrasmall superconducting grains. In that work, 
a modification of the algorithm - called the particle-hole DMRG (phDMRG) - was 
proposed in order to optimally account for the crucial correlations close to the Fermi 
energy. The method was shown to work extremely well for superconducting grains up 
to very large dimensions. The phDMRG was then introduced for possible application 
in large-scale nuclear shell-model calculations |2T> ] 12*7 ] 12%]. 

Several other important applications of the extended methodology have also been 
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reported recently, including applications to electrons in the lowest three Landau levels 

[2H EI3 EB E21 E3] 

The DMRG methodology, both in real space and in the versions appropriate to 
finite Fermi systems, usually violates some symmetries of the underlying Hamiltonian. 
While this is not critical when studying certain features of the associated quantum 
systems, it may be critical when dealing with others. There have been some recent 
efforts to modify the DMRG algorithm jHU E3 EH EZj to permit the preservation of 
symmetries throughout the iterative truncation process. 

In this review, we summarize the rapid progress now being made on the DMRG 
methodology and the wide variety of important applications that have been reported 
recently for strongly correlated finite Fermi systems. We begin in Section II with an 
historical overview of the method and review its application in real space to quantum 
lattices. In the same section, we describe the basic features of the kDMRG extension 
and discuss its application to quantum lattices in more than one dimension. In Section 
III, we begin our focus on finite Fermi systems, discussing the modifications to the 
general DMRG algorithm needed in such cases and reviewing the applications to date. 
Then in Section IV we discuss the issue of symmetry restoration, summarizing recent 
developments and describing in some detail how this new machinery can be used in 
the context of angular momentum restoration. In Section V, we close the presentation 
by providing some perspectives for future applications of the DMRG approach in finite 
Fermi systems. 

2. Historical overview of the DMRG method 

2.1. The many-body quantum problem 

A quantum many-body problem is completely defined by giving a basis of one-body 
states, a Hamiltonian operator and the number of particles. Usually, the one-body 
states are obtained from a self-consistent mean-field calculation, such as Hartree-Fock 
(HF), or by some optimized orbital basis. In atomic and molecular physics 38J, for 
example, the electrons move in the attractive Coulomb field generated by the atomic 
nucleus. The residual electron-electron interaction can be taken into account in mean 
field to give a set single-electron energy states, which then can serve as the starting 
basis. In nuclear physics [39 , the nucleon-nucleon mean field generates a set of self- 
consistent single-nucleon levels, which can provide the basis. Alternatively, we can 
ignore detailed self-consistency and choose a basis for convenience. Convenient single- 
nucleon bases used frequently in nuclear physics include the eigenstates of a harmonic 
oscillator potential (all localized) and the eigenstates of a Woods-Saxon potential (both 
localized and continuum states). 

Though the main goal of this review is to present the latest developments using the 
DMRG method in such finite Fermi systems, we will begin with an historical overview, 
starting with the development of the method for quantum lattices and its applications in 
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that domain. Here, the single-particle basis state is replaced by a single active electron 
level, located at a given site in real space. 

One of the most well-known models of strongly interacting electrons on a lattice is 
the Hubbard model |4U] . It describes the hopping of electrons between nearest neighbor 
sites (atoms) subject to a repulsive Coulomb interaction between them when they occupy 
the same site. The Hamiltonian of the Hubbard model takes the form 



The operator c] a (c ia ) creates (annihilates) an electron on site i with spin projection o. 
The first term represents the hopping of electrons with intensity t. The strength of the 
on-site Coulomb repulsion is U and n ia = c\ a c ia denotes the number operator for an 
electron on site i with spin projection a. 

Assuming that we have a system with N sites, the exact solution of the problem 
amounts to constructing the many-body basis of states in the appropriate Hilbert space, 
calculating the matrix elements of the Hamiltonian in this basis and then diagonalizing. 
In the Hubbard model, each site admits precisely four states, as it can be empty, it can 
contain an electron with spin up, it can contain an electron with spin down, or it can be 
fully occupied by two electrons. Thus, the dimension of the Fock space is naively 4^, 
growing exponentially with the number of sites. Because of conservation of symmetries 
in the many-body basis, e.g. total particle number, total spin, total momentum, etc., the 
Hamiltonian matrix separates into blocks, each with reduced dimension. Nevertheless, 
the order of magnitude of each block dimension is still comparable to the full Fock space 
dimension and thus likewise grows exponentially with the number of lattice sites. For 
this reason, exact large-scale diagonalization using the Lanczos or Davidson algorithms 
is only feasible for a small number of sites. The largest reported exact diagonalization 
of the Hubbard model, for example, was for a 4x4 2D lattice. 

Exact diagonalization, combined with finite size extrapolation, has proven very 
useful when trying to understand the gross features of several lattice models. 
Nevertheless, we would like to do better, by being able to treat larger systems. 
Obviously, if we want to treat larger systems we have to implement an efficient truncation 
strategy, one that permits an accurate description of the low-energy properties of the 
system while only keeping a reasonably small part of the full Hilbert space. 

In the finite Fermi systems that arise in Quantum Chemistry (QC) and in Nuclear 
Structure (NS), similar issues arise as regards the size of the Hilbert space and the 
need for an efficient truncation strategy. In those fields, a variety of methods have 
been developed and used with great success in the study of small systems, such as 
simple molecules and light nuclei. In QC, for example, the Configuration Interaction 
(CI) method [HE] working in a HF basis with a limited set of particle-hole states can 
be used to build a truncated many-body Hilbert space in which to diagonalize the 
Hamiltonian and obtain approximate solutions of the Schrodinger equation. In NS, 
the Shell Model (SM) j^j starts with a nuclear mean field generated by the average 
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interaction of one nucleon with all the others. For simplicity, the mean field is often 
modelled by a harmonic oscillator potential with a spin-orbit coupling, defining a set of 
shells well separated in energy. The nuclear interaction expressed in this basis is typically 
diagonalized for a given nucleus within one major shell for each type of nucleon. Despite 
the enormous truncation inherent in both the CI approach in QC and the SM approach 
in NS, the methods are limited as to the systems that can be treated accurately. In NS, 
for example, exact shell-model diagonalization techniques cannot be applied to nuclei 
beyond the If — 2p shell, leaving outside the scope of detailed shell-model microscopy 
most medium- mass and heavy nuclei. 

Going beyond the limits of exact diagonalization techniques, even within the 
restricted bases just described, requires a new kind of truncation strategy. One that has 
been used with great success in many systems is the Monte Carlo approach. Another 
is the renormalization group (RG) approach |1J, the subject of this review. In this 
method, a set of renormalization transformations that gradually grow the size of the 
system while keeping the dimension of the configuration space constant is introduced. 
The procedure is implemented iteratively. At each step, the system is grown by the 
addition of new degrees of freedom. The enlarged system is then truncated to the same 
size as before, and a renormalization is implemented which transforms the Hamiltonian 
and all other relevant operators to the truncated space. The numerical renormalization 
group of Wilson represented a breakthrough in the application of RG ideas to a non- 
trivial strongly correlated system, the single-impurity Kondo problem Later on 
it became clear that the Wilson RG (WRG) could not be applied to other strongly 
correlated lattice problems [12] . for reasons that we will discuss below. Nevertheless, 
the method and its failure provide the basis for the development of the DMRG. 

2.2. The Wilson Renormalization Group 

In this subsection, we briefly describe the numerical renormalization group of Wilson 
for a generic one- dimensional lattice. While we will have in mind for concreteness 
the Hubbard Hamiltonian the formalism we will describe is general enough to 
incorporate other one- dimensional lattice models, including the Kondo problem. In 
fact, Wilson first formulated his numerical RG after mapping the Kondo problem onto 
a ID lattice in energy space. 

The starting point of a WRG calculation (and as we will discuss later the DMRG 
as well) is a small system for which we know all of the states and for which we can 
calculate the matrix elements of all relevant operators including the Hamiltonian. We 
begin, therefore, by considering a set of L contiguous lattice sites - which we will call 
a block - with L small enough so that we can exactly diagonalize the Hamiltonian in 
the associated subspace. Figure [T] illustrates a chain of lattice sites, with the filled dots 
representing the sites of the initial system - L = 2 in this case - and the open dots 
representing all other sites in the problem. For L = 2, the number of states in the block 
for the Hubbard model is m = 16. 
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Figure 1. Schematic illustration of the sites on a one-dimensional lattice and the 
Wilson Renormalization Group method of enlarging blocks of sites. The sites denoted 
by a bullet • have already been treated, and the next site to the right is in the process 
of being added. 



The WRG procedure begins with a calculation of the matrix elements of all 
relevant operators in the initial block. The relevant operators are those that enter 
the Hamiltonian, namely cj CT , C; CT , rii a , and the Hamiltonian (0) itself. 

Then we continue as follows: 

1. Add the site L + 1 and calculate the matrix elements of all operators in the 
enlarged block space of 4m states. [Note: If we wish, we can partition the states of a 
site and add only 2 states, rather than 4 in each step. When this is done, the enlarged 
block space has 2m states instead of 4m.] 

2. Diagonalize the Hamiltonian and select the lowest m eigenvalues and the 
corresponding eigenvectors. 

3. Transform all operators in the enlarged block to a truncated basis involving these 
m states only. 

4. If L + 1 is lower that the total number of sites N, substitute L + 1 for L and 
repeat steps 1^3. 

The key ingredient in the Wilson renormalization strategy is step 2, whereby m 
states are selected from the 4m (or 2m) states of the enlarged block by a purely energetic 
criterion. Clearly, the accuracy of the method relies on the assumption that the low- 
energy properties of the system are fully contained in the low-energy states of the 
individual blocks into which the system was partitioned. Indeed, this was the case 
for the Kondo problem, which was solved by Wilson by transforming it to a discretized 
energy basis, and for other impurity models like the Anderson model |43.. 

The Kondo model describes the antiferromagnetic interaction of the conduction 
electrons with a localized impurity. Its Hamiltonian in a discretized energy basis takes 
the form 

-1 oo 

H = - (l + A- 1 ) £ A-™/ 2 (f^f n+1 » + fU^U) + 2 Jp/oV^/o, , (2) 

^ n=0 

where /t (fnn) creates (annihilates) an electron in the energy interval n with spin /i, 
J is the strength of the antiferromagnetic interaction between the conduction electrons 
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Figure 2. Schematic illustration of the two-block structure of the Density Matrix 
Renormalization Group. 



and the impurity, p is the density of states and A is a discretization parameter defining 
the intervals of energy. For A > 1 and beginning from the origin (n = 0) of the linear 
chain, the interaction along the chain decays exponentially, assuring the validity of the 
truncation procedure. The WRG was the first numerical implementation of the RG to 
a non-perturbative problem like the Kondo model, for which it had enormous success. 
As such, it is still used to this day to study problems with one or few impurities. 

Unfortunately, the Wilson RG cannot be applied to other lattice problems. For 
example, when it is applied to ID Hubbard models it begins to deviate significantly 
from the exact results in the fourth iteration. Obviously, the problem resides in the fact 
that the truncation strategy was based solely on energy arguments. The solution to this 
problem was proposed by White with his introduction of the DMRG, which we turn to 
in the next subsection. 

2.3. The Density Matrix Truncation Strategy 

The essential limitation of the WRG, as discussed in the previous subsection, is the 
method used for truncation at each RG iteration step. Once the enlarged block is 
constructed by adding the site L + 1 in step 1, the Hamiltonian in the L + 1 - site 
subspace of the total Hilbert space is diagonalized, without considering the coupling 
to the rest of the space. Obviously, the method will work well when this coupling is 
perturbatively small. This is the case for Hamiltonians that can be transformed to the 
form (J2J). In general, however, the WRG performs very poorly for strongly correlated 
lattice systems. In searching for a more suitable method of truncation, White analyzed 
the case of a block immersed in a medium representing the rest of the system, as depicted 
in figure EJ 

Assume that we know the ground state \*f?) of the system and the complete set of p 
states of the block B and t states of the medium M. Then the ground state of the full 
system can be expanded in the product space of the block and the medium as 

l*)=EE *ii\i) B \j) M ■ (3) 

Our aim is to determine the best approximation |^') to the ground state |^) when 
we select only m states of the block, with m < p. To show how this is done, we will 
make use of the concept of the singular value decomposition (SVD) of a matrix jUj. 
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The singular value decomposition of the wave function matrix ^ in (J3J) permits 
its factorization as 

^ ] Q'otU'iaVcij i (4) 
a=X,p 

where the a a are real positive numbers greater than or equal to 0, u is an orthogonal 
square matrix of dimension p and v is a column orthogonal rectangular matrix of 
dimension p x t. The coefficients known as the weight factors of the SVD. 

Up to this point, everything is exact. Now consider what happens if some of the 
weight factors in (@J) are zero or very small. If that is the case, then \1/ can be very well 
approximated by discarding them, i.e. by truncating the sum in (j3J) to those with the 
largest weights. If fact, it can be shown that the approximate wave function factorized 
as in (0J) but retaining the largest m weight factors has the largest overlap with the 
exact wave function for any factorization involving m terms. 

In what follows, we will use the singular value decomposition as a convenient 
factorization for deriving the conditions for optimal truncation in a density matrix 
formalism. It is interesting to note, however, that the same approach has been used 
to develop a powerful new technique for numerically diagonalizing the nuclear shell- 
model problem without actually introducing the DMRG [13] . 

Now consider the reduced density matrix of the block associated with the ground 
state wave function Q, which is defined as 

Pfi = E , (5) 

k=i,t 

with the summation running over the t states of the medium. Inserting the factorization 
(jlj into (|HJ) and using the orthogonality relations of the matrices u and v and the 
normalization of the wave function, we can express the reduced density matrix (jSj) as 

pfj = E a l U iaUaj ■ (6) 

a=l,p 

From this we see that u is the matrix of eigenvectors of the reduced density 
matrix and u a = a 2 a are the corresponding eigenvalues. Each eigenvalue ui a represents 
the probability of finding the block state u a in the ground state of the coupled 
system. The normalization of the wave function \1/ guarantees that the trace of p B is 1. 
Clearly a truncation to m states implies the selection of the largest m eigenvalues and 
corresponding eigenvectors. The discarded p — m states contribute with a probability 
1 — P m = 1 — J2a=m+i u a- This parameter provides an important control on the accuracy 
of the truncation. 

Up to now, we have considered the system as being in a pure state. A similar 
analysis can be carried out for systems in a mixed state. Mixed states naturally 
arise when considering a system at finite temperature (T). They can even arise at 
T = if we decide to target several states with the DMRG truncation strategy. By 
using mixed density matrices associated with several states, we can build a truncation 
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algorithm whereby the gradually enlarging block optimally reflects the physics of all of 
the associated states and not just the physics of the ground state. The most important 
states of the block, when it is immersed in a mixed-state system, are those corresponding 
to the largest eigenvalues of the mixed density matrix, 

pfj = T,w a E n^ k , (7) 

a k=l,t 

where W a is the probability with which a given pure state |\l/ a ) enters in the mixed 
state. In the case of systems at finite temperature, W a is the Boltzmann factor. For a 
mixed state at T=0, we have the freedom to select the weights so as to assign a relative 
importance to each state. 

2-4- The DMRG Method for Quantum Lattices: Real space 

In the previous subsections we described how we can choose the optimal states of a part 
of a quantum system from its density matrix in the presence of a medium approximating 
the rest of the system. Here, we show how this can be used in a systematic treatment 
of many-site quantum lattices. 

We begin by verbally detailing the key ingredients of the method. The basic idea is 
to systematically take into account the physics of all sites on the lattice. The procedure 
begins by starting at the left (or right) edge of the lattice and then gradually growing 
the size of the space that has been taking into account by adding the other lattice 
sites. At each step of this iterative procedure, a truncation is implemented, so as to 
optimally take into account the effect of the most important states within that part 
of the problem. The calculation is carried out as a function of the number of states 
that are maintained after each iteration. This parameter, which we call m, is gradually 
increased and the results are plotted against it. Experience suggests that the results 
converge exponentially with m. Thus, once we achieve changes with increasing m that 
are acceptably small we terminate the calculation. 

In the following, we show how this strategy is implemented at two different levels 
of complexity, referred to as the infinite algorithm and the finite algorithm, respectively. 
We then describe some technical details of importance to the computational feasibility 
of applying the method. 

2. 4-1 ■ The Infinite Algorithm. The infinite algorithm, as we will soon see, can only be 
applied meaningfully in the presence of a uniform (or translationally invariant) lattice. 
Thus, in what follows we assume that the lattice under discussion satisfies this criterion. 

Consider such a lattice, as illustrated in figure EK- Assume that each site on the 
lattice admits I states. In the case of a spin- 1/2 chain, 1 = 2, corresponding to spin up 
or spin down. In the case of the t-J model, / = 3, as the site can also be empty. Finally, 
as noted earlier, in the Hubbard model a site can also have both spin-up and spin-down 
particles and I = 4. 
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Figure 3. Schematic illustration of the infinite DMRG algorithm, (a) The sites on 
a uniform one-dimensional lattice that will be used in subsequent discussion of the 
DMRG algorithm, (b) The two-block structure of the algorithm. A bullet (•) refers to 
a block of sites and a circle (o) to a single site. B refers to the enlarged block, which 
is comprised of bl and si. M refers to the enlarged medium, which is comprised of s2 
and 62. 



We now begin to add sites, by assumption from the left. Until we have included 
enough sites k so that l k > m, no truncation is required. We thus begin our discussion 
after we have added a sufficient number of sites to be in the truncation regime. We 
assume that we have already treated the leftmost k lattice sites and have identified the 
most important m states. 

At that point, we can view the problem as illustrated in figure^. The collection of 
k sites that have already been treated is represented by a • and denoted 61. Immediately 
to its right is another site, which contains I states and is denoted si. We then have 
another site to its right, which also contains I states and is denoted s2. After that we 
have another block, representing the optimal states from sites k + 3 — » 2k + 2. Because 
of the translational invariance of a uniform quantum lattice, all parts of the lattice look 
the same. As a consequence, once we know the optimal states from sites 1 — > k we 
also know the optimal states from any other group of k contiguous sites. This second 
block is denoted 62. And then finally we have other sites to the right that are still to 
be considered. 

We then consider as an enlarged block the m x I states obtained by adding the 
states of site si to those of block 61. We will use the Density Matrix strategy described 
in the previous subsection to truncate from these m x / states to the optimal m states, 
the same number that we had before we enlarged the block. As the medium in which 
to implement the truncation, we use the set of m x I states of the site s2 coupled to the 
block 62. The rationale for growing both the block and the medium at the same time is 
discussed in ref. 
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Following implementation of the density matrix truncation, we make the following 
replacements: 

• The set of m states of the enlarged block now becomes block 61. 

• The next available site to its right becomes si 

• The next site after that becomes s2. 

• The truncated set of m states from the next k + 1 sites becomes block 62. 

This leaves us in the same position as in figureEb, with a block of m states, followed 
by a site of / states, then another site of / states and finally a block of m states. As before, 
we build enlarged blocks out of the first two and out of the second two, respectively. We 
then truncate the enlarged block from 61 and si to the same m states as before. And we 
then use the translational symmetry of the problem to form another set of four objects, 
61, si, s2, 62, and continue the iteration. This iterative process can be continued until 
the block includes half of the sites. At this point, the superblock contains the physics 
of all of the sites and thus gives a representation of the entire problem. 

Such a procedure, in which we grow the system from left to right (or from right to 
left) until all sites have been treated is called the infinite DMRG algorithm. Clearly it 
requires that the lattice be uniform, so that we can use the translational invariance of 
the problem to relate the collective structure of the medium to that of the block. In 
such cases, the accuracy of the procedure still depends on how rapidly the correlations 
between sites fall off. Exact numerical results can be obtained when the correlations are 
predominantly between fairly close neighbors or, more precisely, when the correlation 
length is finite. A typical example is the S = 1 Heisenberg chain 

2.4-2. The Finite Algorithm and Sweeping. Often the above criteria for success of the 
infinite algorithm are not satisfied and it is necessary to do better. When the lattice is 
not unform, there is no natural way to build a medium that properly reflects correlations 
in that part of the lattice. When the correlations do not fall off significantly rapidly, 
there is no way to take into account the effects of distant sites - those that are not part 
of the medium - in defining the collective structure of the block. A way to overcome 
these limitations is to use the so-called finite algorithm, which we now describe. 

The basic idea of the finite algorithm is that we sweep forwards and backwards 
through the set of sites, at each stage updating the information from the previous 
sweep. We will continue to focus here on quantum lattices, and in that context show 
how the method is implemented. We will first carry out the discussion for a uniform 
lattice and then mention the modifications needed when the lattice is not uniform. 

As before we begin by growing the system from the left. For a unform lattice, as 
depicted in fig. EK, we follow the infinite algorithm until half of the sites have been 
treated in the left block and half in the right block. In each iteration we store on hard 
disk all relevant matrix elements within the block for later use. From this point on, 
the superblocks will involve all the N sites of the finite lattice. We then continue to 
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Figure 4. Schematic illustration of the finite DMRG algorithm. B denoted the initial 
block and B' the enlarged block; M denotes the initial medium and M' the enlarged 
medium. 



grow the sites in the left block, one at a time. For each enlarged block, we use for the 
medium the remaining set of sites from the right half. All information on that set of 
sites is obtained by a reflection from the stored information of the left block from an 
earlier iteration, assuming inversion symmetry around the center of the lattice. 

We continue this process until the left block has been grown to include all sites 
but three. At this point, we have stored information on all subgroups of sites, from 2 
through N — 3. That information includes the matrix elements of all sub-operators of 
the Hamiltonian, within the space of states being maintained. 

If the lattice is not uniform, we must use a slightly different strategy, since there is 
no simple prescription for relating the states of the medium to those of the block. One 
possibility is to use the Wilson RG procedure, growing the system one site at a time 
and truncating by energy considerations [19|. Alternatively, one can use the infinite 
algorithm, beginning with blocks at the two ends and gradually expanding them by 
adding sites in between in the direction of the lattice center In either case, we 

continue until all sites have been treated and information on all subgroups has been 
calculated and stored. 

At this point, whether dealing with uniform or non- uniform lattices, we reverse the 
process and sweep from the right, as illustrated in figure El We begin to enlarge the 
right block, by adding one site at a time. The medium for this enlarged block is the 
set of sites on the left, which consists of the block M and the extra site. Note that all 
relevant information on the block M was stored in the previous sweep. 

At this point, we reverse again, sweeping from left to right. The sweep process is 
continued iteratively until the results from one sweep are acceptably close to those from 
the previous sweep. The comparison can be done for the ground-state energy and for 
the density matrix eigenvalues of the retained states. 

Note that the first (or warm-up) stage, in which we sweep through the sites for 
the first time, need not necessarily produce very good results, since we still have the 
opportunity to improve on them in subsequent iterations. It is for this reason that the 
finite algorithm can provide a good approximation to the physics of non-unform lattices. 
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B 



Figure 5. Schematic illustration of block enlargement in the DMRG. A denotes the 
initial block and B denotes the block (or site) that is going to be added to it. 



2.4-3. Some technical details. For the purposes of this discussion, we assume that 
the Hamiltonian describing the system under investigation contains one- and two-body 
terms only and can be written as 

m mnop 

The above Hamiltonian is itself built up out of other operators, which we referred 
to earlier and will continue to refer to as the Hamiltonian sub-operators. They are 

plus their hermitean adjoints. 

A key to the computational feasibility of the DMRG algorithms described earlier 
in this subsection is that it can be readily automated, using information from previous 
steps. Here we discuss this point in somewhat more detail. 

As we have seen, the DMRG approach involves a systematic enlargement of blocks, 
followed by a truncation to the same block size as before the enlargement. Consider a 
representative block enlargement as illustrated in figure EJ The original block A with 
dimension oIa was obtained from an earlier step, at which point the matrix elements 
of all sub-operators of the Hamiltonian (see Eq. Q) were calculated and stored. We 
add to it another block B of dimension oIb, which can either be a block obtained in an 
earlier step or a single lattice site. If it is a block obtained from an earlier step, the 
matrix elements of all sub-operators of the Hamiltonian were calculated and stored at 
that time for the corresponding set of states. If it is a single lattice site, the matrix 
elements of the Hamiltonian sub-operators were stored at the outset of the calculation. 
In either case, the matrix elements of all Hamiltonian sub-operators are accessible in 
storage when the enlargement of the block is implemented. 

At this point, we wish to calculate the matrix elements of the Hamiltonian sub- 
operators in the enlarged block, i. e. in a space of product states 

\iJ) = \i)A \J)b, (10) 

where i — 1, d& and j = 1, ds- 

Consider as a specific example the matrix elements of the one-body operator o) m a n 
in this basis. They can be obtained as a sum over terms involving products of matrix 
elements in space A and space B, viz: 
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(iJla^anlkJ) = A (i\a: m a n \k) A Sjj 
+ B{jW m a n \l)B S i:k 
+ (-) n ^ A (i\aUk)A b(JK\1)b 
-{-T^ A {i\a n \k) AB {j\al\l) B , 

(11) 

where n^ A is the number of particles in state \k) a- Depending on whether the operator 
indices refer to the sites in block A or to those in block B, a different term in the sum 
applies. For example, if both m and n refer to sites in block A, it is the first term that 
applies. Since all of the matrix elements in Eq. (fTTj) were calculated and stored earlier, 
it is straightforward to compute the matrix elements in the enlarged block. 

While we have focused for specificity on the one-body particle-hole operator a^a™, 
it should be clear that the same remarks apply to all of the operators in Eq. (JHJ). 
Obviously, the number of terms in the sum increases with the number of creation and/or 
annihilation operators that define the composite operator being treated. 

Another key step of the DMRG algorithm is renormalization, whereby all sub- 
operators of the Hamiltonian are transformed to the truncated basis of block states 
emerging from the density matrix diagonalization, viz: 

p B \a>=uj*\a> , (12) 

where 

|a>=]TA7|z> . (13) 

i 

The matrix elements of any Hamiltonian sub-operator O in the 4m-dimensional 
enlarged block space \i > can be transformed to the truncated m-dimensional block 
space of states \a > according to 

O a p = Y,XrX^O ij , (14) 

ij 

and then stored for use in a subsequent step of the algorithm. 
2.5. The DMRG Method for Quantum Lattices: Momentum space 

In the previous two subsections, we reviewed the key steps in implementing the DMRG 
algorithm for quantum lattices. The method described worked directly in the real space 
of lattices, growing the system through the addition of sites. 

A key to the success of the DMRG algorithm is being able to define an appropriate 
order in which to add layers of the onion. In the case of lattices, this means an 
appropriate order in which to add sites. In one dimension, there is a natural order, 
namely to add sites either from the right or from the left, always adding the next 
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available "nearest neighbor" site(s). This is natural since the interactions that couple 
sites tend to fall off rapidly with increasing distance between them. But what about 
lattices in more than one dimension? As an example, consider a 2D square lattice. Near 
any given site there are four "nearest neighbor" sites. What is the correct path through 
this set of lattice sites to consider? 

For such reasons, it may be useful when dealing with lattices in more than one 
dimension to consider an alternative labelling of the states of the system, one that is 
more conducive to a natural ordering. 

In 1994, Xiang [TH] proposed the possibility of carrying out DMRG calculations in 
momentum space, rather than in real space. In this section, we review the new ideas that 
emerged from this proposal. As we will see in the next section, this idea provides the 
inspiration for the algorithms that have since been developed for finite Fermi systems, 
the principal focus of this review. 

To make the ideas most transparent, we will carry out the discussion in terms of 
the two-dimensional Hubbard model with nearest-neighbor hopping, as was done in the 
original paper of Xiang. The Hamiltonian for this model was given earlier in Eq. ([!}. 

The operator that creates a particle in Bloch state kcr with spin a and momentum 
k and the one that creates a particle with the same spin a but at lattice site i are related 
according to 

^L = 4rEe ikri 4, (15) 



'N 

where N is the total number of sites in the lattice and r; is the coordinate-space location 
of site i. 

Using this, it is straightforward to transform the Hubbard Hamiltonian (JJJ to 
momentum space, 

H = e k^k<x + ^ E C p-qT C k+q|Ck|C pT , (16) 

kcr P,k,q 

where 

= cl a c ka (17) 
is the number operator for particles with spin a and momentum k and 

e k = -t e ik < ri - T ^ (18) 

<i,j> 

is the energy dispersion of the particles, which is explicitly diagonal in the momentum 
basis. 

In particular, by diagonalizing the one-body hopping term for a 2D lattice we find 

that 



ek = — 2t (cos k x + cos k y ) , 
k = (2-Kn x /L , 2irn y /L) , 

n x , n y = -L/2 + l,...,L/2 ; L = y/W . (19) 
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Table 1. Energy dispersion for the 2D 4x4 Hubbard model. 



e/t Degeneracy 



(k x ,ky) 



-4 
-2 

2 
4 



1 

4 
G 
4 
1 



(0,0) 

(-rr/2,0), (0,-7r/2), (tt/2,0), (0,tt/2) 
(-tt/WU (-tt/WU (ir/2,-n/2), (n/2,ir/2), (tt,0), (O.tt) 
(-7r/2,7r), (7r/2,7r), (7r,-7r/2), (^vr/2) 

(7T,7r) 



In table ^ we give the momentum "spectrum" for a 2D 4x4 lattice. Possible 
values for k x and k y in this case are -tt/2, 0, 7r/2, and n. What is clear is that there is a 
large degeneracy at half filling. While it is natural to consider filling momentum-space 
levels in order of increasing energy dispersion, it is clear that there will be a problem 
with degeneracies (or near degeneracies when more general hopping terms are used). 
Thus, the rationale we gave above for introducing momentum space - namely to provide 
a more natural ordering of levels (onion layers) - may not turn out to be warranted. 
Nevertheless, we will discuss what happens when the methodology is in fact used. Before 
doing so, however, we will first discuss some issues related to the implementation of the 
momentum-space DMRG algorithm, which we subsequently refer to as the kDMRG. 

2.5.1. The size of the problem. As emphasized by Xiang in his original paper, the 
two-body interaction appearing in (fTT)|) contains of order N 3 terms, since each of the 
sums goes over N possible momentum values. Thus, a straightforward application of the 
DMRG algorithm in momentum space would require iteratively calculating and storing 
matrix elements of 0(N 3 ) operators, which would obviously be very time and memory 
consuming for a large number of sites. Xiang figured out a very clever way to group the 
operators, so as to dramatically reduce the number whose matrix elements need to be 
calculated and stored. Since an analogous strategy will be required in all extensions of 
the usual rDMRG method, we will now discuss how this is accomplished. 

Consider the following set of operators, operating in a subspace A of the full 2D 
Hubbard Hilbert space: 





a iG°> °) = £ 4(ft <r)ao(p + ft o) 



a 2{p) = a o(ft TKO + q, I) 



fl3(p, <*) = X] o!(?i> -o")ao(?2, -o>o0 + ft. - ft;, <r) (p, °) £ A 



a i(p) = Yl a o(ft IHO - ft t) • 



(20) 
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The total number of such composite operators is 6N, significantly less than the total 
number of sub-operators of the full Hubbard Hamiltonian in momentum space. 

For notational purposes, we only use the above set of composite operators when 
referring to block A. When considering another block B, we systematically replace a, 
by h. 

The key point is that all matrix elements needed for the iterative algorithm can be 
obtained from these composite operators. To see this, consider for example the Hubbard 
Hamiltonian in a product space AxB. It can be written as 

H AxB = H A + H B + ^ Yl jUt^ife a ') b i(-p, + 6j(p, a')a 3 (p, a') 

+(a <-►&)] + 4(p)h(p) - b 2 (p)al(p)} + h.c. (21) 

Note that in the A space, we only need the set of a operators and correspondingly in 
the B space we only need the analogous set of b operators. 

In more complex spaces, with larger numbers of coupled blocks, the same remarks 
can be straightforwardly shown to apply. 

The above development is specific to the Hubbard Hamiltonian of Eq. (0). However, 
as we will discuss later, similar considerations are essential to preserve the computational 
feasibility of the DMRG algorithm in the more general quantum problems that arise in 
the description of finite Fermi systems. 

2.5.2. Translational invariance and momentum conservation. As noted earlier, the 
real-space DMRG when applied to uniform quantum lattices benefits greatly from the 
fact that all positions on the lattice look the same. Thus, once a certain portion of the 
lattice has been treated, the same results can be transferred to any equivalent set of 
sites. 

The momentum-space DMRG loses this simplicity. A given group of quantized 
momentum levels has no obvious relation to any other group. 

Of course the equivalence of a given group of sites on a uniform lattice to any other 
analogous group was a consequence of the translational invariance of the problem. That 
symmetry is not lost by transforming to momentum space; it simply shows up elsewhere. 
In particular, in momentum space the translational invariance of the lattice shows up 
in terms of momentum conservation. When building the superblock, we need not mix 
states of different total momentum. This has the advantage of reducing the size of the 
superblocks in momentum space compared to real space, even for the same assumed 
block dimension. 

In building states of a given total momentum in the superblock, we must maintain 
all allowed momentum states for the blocks themselves. The resulting density matrix 
for the block will thus contain different momenta. However, it will be block diagonal 
in momentum and thus likewise be effectively reduced in dimension relative to the 
corresponding matrix in real space. 
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Similar considerations play a role in applications of the DMRG to finite Fermi 
systems, where symmetries likewise must be treated carefully. In rotationally-invariant 
systems, for example, angular momentum is a conserved quantum number. As we will 
discuss in Section IV, full incorporation of angular momentum symmetry can be built 
into the formalism, but at an increase in complexity. Minimally, however, when working 
in a basis of single-particle states with well-defined angular momentum projection, 
we can build states with conserved m-projection and thereby reduce the size of the 
superblocks that must be diagonalized. Any additive quantum number associated with 
an Abelian symmetry can be used in this way to reduce the size of the superblock and 
to make the reduced density matrices block diagonal. 

2.5.3. The warm-up phase of the kDMRG algorithm. In a real-space description of 
uniform lattices, we can use the equivalence of contiguous sites throughout the lattice 
to facilitate the introduction of a medium with which to build a DMRG superblock in 
the warm-up phase. In momentum-space, as noted above, this is no longer the case. On 
the one hand, this makes it clear that it is not appropriate to use the infinite DMRG 
algorithm when working in momentum space. But even with the finite algorithm it 
is necessary to adopt a strategy for implementing the warm-up phase, where the first 
approximation to the description of the the various size blocks is obtained. 

In the work of Xiang and the others who have implemented the kDMRG, the 
procedure has been to perform truncation in the warm-up phase according to Wilson's 
numerical RG procedure, namely to truncate to the lowest eigenstates of the enlarged 
block Hamiltonian. As for a real-space description of non-uniform lattices, the 
expectation is that subsequent sweeps should correct any deficiencies that might exist 
in this procedure. 

2.5.4- Applications of the kDMRG. Having described the key elements of the kDMRG 
method, we now discuss some of the applications that have been reported. 

We first discuss some of the results reported by Xiang in the paper that first 
introduced the method [IHj. The calculations were carried out for the ID and 2D 
Hubbard models, using the standard Hamiltonian of Eq. (0). They were performed 
as a function of the number of lattice sites N and the ratio U/t of the strength of the 
on-site Coulomb repulsion to that of the one-body hopping term. In the case of the ID 
model, all calculations were done at half filling, whereas the 2D calculations considered 
different filling fractions. In all cases, the number of states m kept in each truncation 
was 1000 and comparison with exact results was reported where available. 

Applications in ID showed two important points: (1) that the momentum-space 
DMRG method is significantly more accurate than the Wilson RG when applied to the 
ID Hubbard model, and (2) that the kDMRG results for ID lattices are less precise 
than those for real-space DMRG. 

The latter conclusion is worthy of further comment. As noted earlier, the accuracy 
of the DMRG algorithm depends critically on how rapidly the correlations fall off. In 
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Table 2. Comparison of the ground state energy calculated for the 4x4 Hubbard 
model. The exact results are compared with those obtained using the kDMRG method, 
the Quantum Monte Carlo (QMC) method and the Stochastic Diagonalization (SQ) 
method for different numbers of electrons N and different values of U/t. From Table 
II of ref . | 19| . Reprinted by permission of the American Physical Society. 



U/t 


N 


Exact 


kDMRG 


QMC 


SQ 


2 


16 


-18.01757 


-18.012 






4 


14 


-15.74459 


-15.673 






4 


16 


-13.62185 


-13.571 


-13.6 


-13.59 


8 


16 


-8.263 


-8.48 







real space, the correlations decay with distance since the Hamiltonian only couples 
neighboring lattice sites. An analogous statement does not apply in momentum space, 
where the interaction is highly non-local and can strongly couple momentum states 
that are far apart. A critical question there is how strong is this coupling relative to 
the single-particle splittings between the Bloch states. Nevertheless, the conclusion to 
draw from this work is that when studying ID lattices using the DMRG algorithm, the 
real-space formalism remains the method of choice. 

Xiang's results for the Hubbard model on a 4 x 4 2D lattice are summarized in table 
121 The results are shown for several values of U/t and for several numbers of electrons. 
Where available they are also compared with results of Quantum Monte Carlo (QMC) 
and Stochastic Diagonalization (SQ) calculations. There are two key points to note. 
The first is that the method works very well for weak on-site repulsion. This should 
not be surprising, since it is exact when {7 = 0. In that case, the exact solution - a 
Slater determinant of Bloch states up to the Fermi momentum - is a trivial one-state 
kDMRG wave function. The second point to note is that the method continues to work 
reasonably well even for stronger coupling, producing results that are competitive with 
those obtained using the two other popular many-body methods against which it was 
compared. 

Several years later, Nishimoto and collaborators [H)] reinvestigated the usefulness 
of the kDMRG method for treating the Hubbard model. They redid some of 
the calculations reported by Xiang and also carried out a much more systematic 
investigation, considering several different hopping scenarios. An important difference of 
their calculations relative to those of Xiang was in the blocking scheme. These authors 
applied White's blocking scheme, adding two sites at a time to enlarge the superblock. 
Xiang, in his kDMRG work, added only one site at a time. This makes it somewhat 
difficult to draw comparative conclusions, as the sizes of the superblocks that would arise 
in the two blocking approaches are different for the same number of states m retained 
in each block. 

We summarize here some of the key conclusions of ref. [IE]. 
• In ID with nearest-neighbor hopping, Xiang's results do not match theirs. Several 
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possible explanations were put forth, including that Xiang's calculations may not 
have converged and that it could be a consequence of the different blocking schemes 
used in the two calculations. 

• In 2D, both the rDMRG and the kDMRG become less accurate with increasing 
dimensionality, albeit for different reasons. In real space, it relates to the level of 
off-diagonal coupling. In momentum space, it relates to the increased degeneracies 
and to the long-range nature of the interaction. 

• The kDMRG is less dependent on the form of the hopping term. In real space, the 
method becomes rapidly worse as the hopping range increases; this is not the case 
in momentum space. 

• In agreement with Xiang, they find that the accuracy of the DMRG method 
becomes worse with increasing strength U, scaled by the band width. 

• In the case of the 4x4 2D lattice at half filling, the kDMRG was found to be more 
accurate than the rDMRG for U < 8t. 

3. DMRG for finite Fermi systems 

The real-space DMRG has proven extremely successful in describing the physics of one- 
dimensional systems with short-range interactions. As discussed in Sect. II. B. 4, the 
method can be extended to 2D systems by selecting a particular path covering the lattice 
at the expense of introducing long-range interactions. These long-range interactions 
arise between neighboring sites in the 2D lattice, which may be located far away within 
the selected ID path. The interplay between short-range and long-range interactions 
works to the detriment of the optimal accuracy of the ID DMRG and to an increase 
in the memory requirements to implement the algorithm, because of the necessity of 
storing the matrix elements of long-range correlators. In practice, real-space 2D DMRG 
calculations were therefore restricted to small clusters or to ladders up to six legs 

The kDMRG represented a possible way to extend the DMRG methodology to 
larger dimensions. Unfortunately, the results obtained for the Hubbard models of 2D 
lattices were not very encouraging and the procedure was not pursued very much. 
However, the formal procedure developed by Xiang to reduce the number of matrix 
elements needed to store at each iteration was a critical step in the subsequent 
application of the DMRG methodology to finite Fermi systems, to which we now turn. 

In finite Fermi systems one typically starts with a one- plus two-body Hamiltonian 
expressed in some single-particle basis, 

H = J2 T iA c i + E Vijkic\c]c lCk , (22) 

ij ijkl 

where takes into account the kinetic energy and the mean-field interaction between 
particles and Vijki represents the residual interaction. 

A typical basis used in studies in quantum chemistry, nuclear physics, matter 
waves, etc. involves the use of the Hartree-Fock (HF) mean-field approximation. Other 
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frequently-considered possibilities include the Kohn-Sham basis or approximate mean 
fields. As an example of the latter, nuclear structure calculations often use for simplicity 
a 3D harmonic oscillator basis or a basis comprised of the eigenstates of a Woods-Saxon 
potential, as approximations to the nuclear mean-field potential. 

The first attempt to apply the DMRG to a 3D finite Fermi system was performed 
by White and Martin [20] • They proposed the use of the DMRG for ab initio calculations 
of the electronic structure of molecules. The method was further studied and refined 
for calculations in quantum chemistry and for application to quantum Hall systems 

HUB I!!! Ell EH] ■ Other forms of the DMRG procedure were developed and applied 
to ultrasmall superconducting grains [2U EE] and the nuclear shell model PHI 123 EE]- In 
the following subsections, we describe the progress achieved in each of these domains, 
as well as the difficulties and limitations of the procedures that were used. 

3.1. Quantum Chemistry 

Ab initio calculations in Quantum Chemistry (QC) are based on the approximate 
treatment of a system of iV electrons moving in a field created by M fixed atomic 
nuclei and interacting through the Coulomb force. The electronic Hamiltonian of the 
system is 

N ^2 N,M z ^ e 2 ^ N e 2 

i=\ * i=i, a =l I 1 a ' *tyj=x\'t ' 3\ 

There are two possible avenues |HB] for the treatment of the electronic Hamiltonian 
(J23|) . One involves the use of density functional theory (DFT), which has been used 
extensively and with great success in the treatment of fairly large molecules. The other 
involves the use of Hartree-Fock (HF) mean-field theory as a starting point for an 
improved treatment of the electron-electron interaction. 

In principle, the best means of treating the electron-electron interaction in a 
HF framework is through the full Configuration Interaction (CI) approach |3S]. In 
this approach, the electron-electron interaction is exactly diagonalized in a basis of 
orthogonal states generated by the HF method. For problems that are not too large, 
this can be implemented through the use of diagonalization routines appropriate to 
sparse matrices, such as Lanczos or Davidson. With present computational resources, 
matrices of dimension 10 7 to 10 s can be treated in this way, with little hope of a dramatic 
increase in the future. Other approximate but very accurate methods to treat the 
residual electronic interactions, such as Coupled Cluster Theory (CCT), are often used, 
but like the full CI they scale exponentially with the number of basis states, thereby 
imposing limitations on the Hamiltonian matrix dimensions that can be treated and 
permitting complete calculations for relatively small molecules only. 

The DMRG method was proposed by White and Martin |2*U1 |2*T| as an alternative 
procedure for approximately diagonalizing the residual electronic interaction, but in a 
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much more computationally efficient way. Despite requiring a much slower polynomial 
computational effort, the method nevertheless has been shown to provide highly accurate 
results. More precisely, the DMRG computational effort scales as 0(n 4 m 2 ) where n is 
the total number of basis states and m is the number of states retained in the iterative 
algorithm. This suggests that application of the DMRG algorithm to the field of QC 
could open up the possibility of a wide range of highly accurate molecular structure 
calculations, well beyond the limits of the full CI. 

The usual way to generate a second quantized Hamiltonian in a finite molecular 
basis is to first consider a small number of atomic orbital states for each atomic nucleus, 
generating a local non-orthogonal basis, and then to use HF mean-field theory to produce 
an orthogonal molecular basis of n orbitals, doubly-degenerate in the spin quantum 
number. The second- quantized Hamiltonian in this basis is 

H = J2 T iJ c t c ^ + \ Yl Vijkiclc] a ,ci a/ c ka , (24) 

ija ijklaa' 

where c\ a (ci a ) is the creation (annihilation) operator of an electron in orbital i with 
spin projection a. The matrix T is a radial integral of the electron kinetic energy and 
the electron-nucleus interaction in the HF basis, while the matrix V is a corresponding 
radial integral of the repulsive Coulomb interaction between electrons. 

The key point in applying the DMRG procedure to QC problems is to treat each 
orbital as a "site" in a one- dimensional lattice. This site can be empty, singly occupied 
with spin up, singly occupied with spin down, or doubly occupied. The situation is 
analogous to the kDMRG and in fact it was the use of the Xiang trick for reducing 
the three-index operators to an effective one-index operator that made the procedure 
tractable in QC. The most complex operators that need to be calculated and stored 
in each iteration are those with two electronic creation and/or annihilation operators, 
requiring storage space of order n 2 m 2 . 

As we will see below, a key issue, which is still not in our view fully resolved, is 
how to establish the optimal relative order among the orbitals. Two possibilities are 
to order the orbitals in terms of HF energies or to maximize the interactions between 
neighboring orbitals. A satisfactory resolution of this problem will ultimately lead to 
an optimized DMRG protocol for applications in QC. 

As an illustration of the first tests of the method, we show in fig. El DMRG results 
for the methane molecule (CH4) [21] ■ The system consists of 10 electrons moving in a 
minimal basis of 9 orbitals. The one- and two-body matrices in the Hamiltonian (|24j) 
were obtained following a HF calculation. The HF molecular orbits played the role 
of sites in a finite one-dimensional lattice and the finite DMRG algorithm was applied 
to the problem, with the ground state (GS) targeted. The DMRG results for the GS 
energy are plotted as a function of the number of states m retained in the iterative 
procedure. For comparison, the figure also includes the CAS(/i, v) results corresponding 
to diagonalization of the Hamiltonian ()24j) in a Complete Active Space of v valence 
orbitals and \i active electrons. 
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Figure 6. Ground-state energy E of CH4 as a function of the number of states retained 
m at each DMRG iteration. The lines are the Hartree-Fock and CAS results. Taken 
from [2]. Reprinted by permission of John Wiley & Sons, Inc. 



The DMRG results exponentially converge to the exact results of CAS(10, 9) in 
which the complete active space of 9 molecular orbits is taken into account. For rather 
small m values (~ 50) the DMRG results are already comparable to CAS (8, 8) in which 
two electrons are frozen in the lowest HF orbital. 

The influence of the choice of basis on the global procedure was also investigated 
in (21], with a focus on the singlet state of the H-He molecule. The standard definition 
of the molecular orbits based on HF mean field approximation was compared with the 
DFT definition of natural orbitals in two different approximations - the Local Density 
Approximation and the Generalized Gradient Approximation. The results obtained 
were fairly insensitive to the choice of basis, showing similar convergence properties. 

The efficiency of the DMRG procedure in QC has been further studied in ref. 
|22j . Though the method was found to perform very well for the molecules Be2 and 
HF, comparison with full CI results for the N2 molecule did not produce satisfactory 
results. The reasons for this failure were not conclusive, but there were indications that 
a non-optimal ordering of orbitals may have played an important role. 

A very careful study of all of the ingredients that define the performance and 
efficiency of the DMRG algorithm in QC was reported in ref. [22]; where the reader 
can also find a detailed and pedagogical presentation of the method. The subject of the 
orbital ordering was addressed in this work, with the idea of minimizing the range of 
the interaction matrix elements in order to render the correlation length of the system 
a minimum. After trying various prescriptions, the authors came to the conclusion that 
the simplest method was to minimize the band-with of the one-body matrix elements 



The Density Matrix Renormalization Group for finite Fermi systems 



25 



Table 3. DMRG results for the water molecule for various values of to, compared 
with three levels of approximation in CCT. Taken from |23| . Reprinted by permission 
of the American Institute of Physics. 



m 


E/H 


SE/mH 


100 


-76.2545 


2.1 


200 


-76.2559 


0.71 


300 


-76.25632 


0.32 


400 


-76.256477 


0.157 


500 


-76.256540 


0.094 


600 


-76.256592 


0.042 


750 


-76.256617 


0.017 


900 


-76.256624 


0.004 


CCT 


-76.252503 


4.131 


CCT(T) 


-76.255907 


0.727 


CCT(TQ) 


-76.256846 


0.202 



Tij. They found that the symmetric reverse Cuthill-McKee reordering method could 
bring the T matrix to a block diagonal form, avoiding long-range hopping terms. 

The example of the water molecule, which had been treated by White and Martin 
in [20] as a fi rs t test of the DMRG method in QC, was re- investigated in [2H], to compare 
the accuracy of the DMRG not only with a full CI calculation but also with current CCT 
approximations. This problem involves 8 electrons moving in 25 active orbits (the Is 
orbit was frozen). We reproduce in table HH their results. The exact GS energy given by a 
complete CI diagonalization is -76.256634#. CCT(T) and CCT(TQ) refer to Coupled 
Cluster calculations with triples and with triples and quadruples, respectively. It can be 
seen from the results that for m = 400 the DMRG already has comparable accuracy to 
CCT(TQ). Furthermore, and most importantly, since DMRG is a variational theory we 
can be sure that the numerical result it yields is an upper bound to the GS energy. This 
is not the case for CCT, as seen from the CCT(TQ) result for which there is overbinding 
of 0.2^H. 

Analysis of larger molecules, such as C2H4, indicates that it is difficult for the 
DMRG method to capture long-range electron correlations because of the sequential 
procedure it uses for adding orbitals. Nevertheless, the fact that the DMRG results 
obtained in ref. were comparable to the CCT results including quadruples suggests 
that it should be taken seriously in such calculations, with possible improvements now 
to be discussed. 

A difficulty in drawing firm conclusions from the studies just discussed is that exact 
results from full CI calculations were not available in all cases, making it impossible 
to confidently assess the true accuracy of the DMRG and CCT calculations. One of 
the control parameters of the DMRG method is the number of states m retained in 
the algorithm. The variational properties of the DMRG assure that the GS energy 
decreases, or at worst stays constant, with increasing m values. Indeed, in many 
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Figure 7. Dynamically selected number of left and right blocks states for two values of 
the minimum threshold m min — 16, 64 and the relative error as a function of iteration 
for three values of m m in = 16, 64, 164. In all cases the 10 electrons of the H2O molecule 
were correlated in the double-zeta water model with 14 orbitals, and TRE max = 10 -10 
was set in advance of the calculations. Taken from ref. Reprinted by permission 
of the American Physical Society. 

cases an exponential or modified exponential j2H| convergence with m has been found. 
Nevertheless, even in such cases it is difficult to correlate the number of states m with the 
final precision of the procedure. In ref. jSI], a new method for determining the accuracy 
of the DMRG in advance, the Dynamical Block State Selection (DBSS) method, was 
proposed. The truncation error (TRE) in each iteration is TRE = 1 — J2a u a, where 
u a are the eigenvalues of the reduced density of the block in descending order. Keeping 
the same number m of states in each iteration induces fluctuations in the TRE, making 
it impossible to infer the final accuracy of the procedure. The DBSS method inverts 
the procedure for selecting states by defining a maximum truncation error TRE max at 
the start, which amounts to a dynamical selection of the number of states kept in each 
iteration. The DBSS is supplemented by the selection of a lower threshold m m i n for the 
number of states retained. The first and second panels of figure extracted from ref. 
|51j . show the number of states dynamically selected in each iteration for two values of 
m m i n in a calculation for a water molecule with 10 electrons and 14 orbitals. The third 
panel shows the relative error in the GS energy, defined as (E DMRG — Ea)/ Ed- 
it is interesting to note that even though the maximum number of states does not 
depend on m min the rate of convergence is faster for larger m m i n values. In the example 
of the water molecule, the HF single particle levels were ordered with increasing values 
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of energy. In other test cases, e.g. for the CH 2 molecule j^I], this ordering led to 
a local minimum and only by changing the ordering was it possible to obtain global 
convergence. 

The problem of the single-particle ordering was investigated recently in ref. |52j . 
where a protocol for the initialization of the procedure and for an optimal ordering of 
the levels was proposed. Using concepts borrowed from Quantum Information Theory, 
it was shown that the optimal ordering corresponds to locating the levels with maximum 
entropy at the center of the chain. This ordering is not always equivalent to the Cuthill- 
McKee ordering that had been proposed in |23j . 

There have been several other applications of the DMRG algorithm in Quantum 
Chemistry to large polymers jS3 EH E3 E3 EH EH EH1 EH EI] • These applications have 
typically utilized the original real-space DMRG algorithm building on the fact that 
polymers can be represented either by a chain or by another cyclic structure. Since our 
focus is on applications that require the machinery appropriate to finite Fermi systems, 
we will not discuss these applications here. 

3.2. Ultrasmall superconducting grains 

In 1959, Anderson jH2j posed the question of what would be the lower size limit for 
a metallic grain to sustain superconductivity. Since the average level spacing d in 
a metallic particle is inversely proportional to its volume, for small enough particles 
(namely when d is of the order of the bulk gap A) the scattering of pairs across the 
Fermi surface ceases to be an efficient mechanism for lowering the energy of the system. 
At this point superconductivity should therefore disappear. Anderson estimated the 
critical grain size to be in the range of nanometers. This old question has recently 
been revived in the context of experiments performed by Ralph, Black and Tinkham 
inB] on ultrasmall Aluminum grains. They succeeded in constructing a Single Electron 
Transistor whose central island was an Al grain of a few nanometers with an estimated 
mean level spacing d ~ OAmev, of the order of the bulk gap. The experiments showed 
the existence of a spectroscopic gap that could be driven to zero by the application of 
magnetic fields. 

The Hamiltonian proposed in [HI] to study the physics of such small grains is of 
the reduced BCS form, 

q n 

ff = EE( £ 3- ^) C )<r C 3° ~ Xd J2 4+4-Cj-C j+ , (25) 
j=l cr=± i,j=l 

where i,j = 1,2, label the single-particle energy levels, which are assumed to 

be equally spaced (ej = jd), and the electron annihilation operators associated 

with the two time-reversed states \j,a = +) and \j,o~ = —). Finally, \i is the chemical 
potential and A is an adimensional BCS coupling constant, whose appropriate value for 
Al grains is 0.224 [S3]. The systems described by this Hamiltonian have been treated 
at half filling for which the number of pairs iV equals the number of levels Q. 
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It soon became clear that the standard BCS approximation, which so accurately 
describes superconductivity for sufficiently large systems, is not well suited to treat 
the systems of physical interest, for which iV ~ 100. In such systems, there are 
strong number fluctuations proportional to y/N. Several levels of approximations were 
thus applied to the Hamiltonian (|23|) . to see whether any could describe the physics 
of pairing in finite systems. In increasing order of complexity, they included i) BCS 
approximation with parity projection [S3], ii) Lanczos diagonalization up to Q = 23 
with finite-size corrections [66j, and iii) number-projected BCS jUZl approximation. All 
of these approximations predicted a phase transition from superconducting behavior 
to a normal state at a critical grain size, and in all cases the transition was abrupt. 
In contrast, the DMRG method - when applied in a somewhat modified form [23] 
showed that the crossover as a function of decreasing grain size had a smooth reduction 
in pairing correlations. 

Following completion of [23j, the authors learned that the Hamiltonian ()25|) had 
been solved exactly by Richardson during the 1960s [SHJ - This made it possible to 
carry out a comparison of the DMRG results with those obtained from Richardson's 
exact solution. For the largest system treated, with 400 doubly-degenerate levels, the 
condensation energy was E% MRG = —22.5168 with a predicted relative error of 10 -4 . 
Richardson's exact result was E^ xact = —22.5183141, confirming the accuracy of the 
phDMRG. and showing that the relative error predicted by the DMRG method was 
in fact appropriate. On this basis, it is clear that the DMRG method provides an 
essentially exact description of the physics of pairing in finite systems and that the 
correct crossover from superconducting to the normal state is indeed smooth and not 
abrupt. 

The modifications implemented in the DMRG method for its application to small 
metallic grains were based on the fact that in a finite Fermi system the Fermi energy 
(Ep) naturally separates the single-particle space into a primarily occupied subspace and 
another primarily empty subspace. The method that emerges from such a construction 
was called the particle- hole DMRG (phDMRG), which we now briefly describe. 

Consider the set of single-particle levels shown in fig. [HJ Note that they are divided 
into two sets. The particle levels are those above the Fermi surface, which should be 
essentially empty. The hole levels are those below the Fermi surface and they should be 
essentially filled. In the phDMRG method, one first defines particle and hole blocks from 
the levels nearest to the Fermi surface, and then gradually grows those blocks by adding 
levels progressively further away. The figure is meant to signify that we have already 
created particle and hole blocks involving the first two available levels, respectively, and 
are enlarging them by adding the third level(s). The particle block then serves as the 
medium for the hole block and the hole block as the medium for the particle block as we 
implement the DMRG truncation strategy. This process is continued until all particle 
and hole levels have been treated. 

The phDMRG procedure, as just described, is related to the infinite DMRG 
algorithm in that sweeping is not implemented. It is extremely difficult to implement 
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particle levels 



hole levels 



Figure 8. Schematic illustration of the phDMRG method for finite Fermi systems. 
The single-particle levels are divided into primarily occupied hole levels and primarily 
empty particle levels. The first two levels above the Fermi energy Ep have been formed 
into a particle block and the first two below have been formed into a hole block. The 
third levels above and below are being added to form enlarged blocks. 



sweeping with particle and hole blocks while at the same time preserving the four block 
structure of the method. As such, it is expected that the phDMRG will only work 
well for systems where the correlations fall off rapidly as we move away from the Fermi 
surface. The key point is that at the beginning of the iterative procedure the levels 
that are outside the particle-hole superblock do not impact on the selection of the 
optimal state (s) with which to truncate the two blocks. This is a recurring limitation 
of the infinite algorithm, which can be partially overcome by combining the DMRG 
truncation strategy with the use of effective interaction theory. In the application to 
ultrasmall superconducting grains, this was implemented through a phenomenological 
renormalization of the pairing strength at each iteration. When the system described 
by the pairing Hamiltonian ()25j) for Q levels at half filling is considered in BCS 
approximation, the bulk gap is given by = G?f2/2sinh(l/A). The effect of the levels 
outside a superblock that consists of 2n levels can thus be taken into account in an 
average way by introducing an effective coupling strength A n - corresponding to step n 
of the iteration - such that the gap in the superblock being considered equals the bulk 
gap. This is guaranteed if the renormalized coupling strength A n is chosen according to 

2(n + 1) 

sinhA n = 1 ^ J smh\ . (26) 

Note that the pairing strength assumes its full value when all levels are included, i.e., 
when n = ^p. 

The effect of renormalizing the coupling strength in this way can be seen in table 
HI where we show results for the condensation energy in units of the mean level spacing, 



The Density Matrix Renormalization Group for finite Fermi systems 



30 



Table 4. The effect of renormalization on the results of phDMRG calculations for 
ultrasmall superconducting grains. 



m 




-^dressed 1 & 


dim'Hioo 


40 


-40.32502 


-40.49884 


1246 


50 


-40.44623 


-40.50014 


2108 


60 


-40.46887 


-40.50061 


3032 


70 


-40.48878 


-40.50068 


3622 


80 


-40.49588 


-40.50072 


4820 


90 


-40.49815 


-40.50074 


6306 


100 


-40.49919 


-40.50075 


7778 


110 


-40.49983 


-40.50075 


9720 



without renormalization E^ are and with renormalization E^ ressed , for different values 
of the number of states m maintained. The system considered has Q = 100 at half 
filling and the bare coupling constant is A = 0.4, for which the system is well into 
the superconducting phase. In the last column, we show the dimension of the largest 
superblock required in the iterative process, which should be compared with the total 
dimension of the problem, 10 29 . 

As can be seen from the table, there is a significant improvement in convergence 
when renormalization of the coupling constant is implemented in the phDMRG 
algorithm. For m = 100, for which the largest matrix diagonalized had a dimension 
of roughly 8, 000, the condensation energy (column 3) has converged to seven significant 
figures, while the standard phDMRG procedure (column 2) without renormalization 
is still far from convergence. The results of the complete calculation for ultrasmall 
superconducting grains are shown in fig. 01 The condensation energies for even and 
odd systems are displayed as a function of the mean level spacing in units of the bulk 
gap. The solid lines correspond to the numerically exact phDMRG calculation and 
the squares and triangles, for even and odd systems, respectively, correspond to the 
PBCS calculations of jSS]- The phDMRG results were obtained with m = 60 with an 
estimated relative error for the largest system considered (Q = 400 levels) of 10 -4 . The 
phDMRG results are consistently below the PBCS results, representing a considerable 
gain in correlation energy. Most importantly, while the PBCS results suggest an abrupt 
change of the condensation energy at a critical value of the level spacing, the phDMRG 
results unveiled the true nature of the transition. There is a very smooth crossover 
as the size of the grains decrease, from superconducting grains to metallic grains with 
strong pairing fluctuations. 

The fact that BCS approximation ceases to provide a good approximation for small 
enough grains suggests that other phenomena associated with superconductivity may 
also be strongly modified in that regime. One example is the Josephson effect, whereby 
correlated Cooper pairs tunnel from one superconducting system to another due to the 
differences in phases of their respective order parameters. When BCS approximation 
breaks down, it is no longer meaningful to speak of superconducting phases, raising the 



The Density M 

-1.5 

-2.0 

-2.5 
< -3.0 
m -3.5 

-4.0 

-4.5 

-5.0 

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 

d/A 

Figure 9. Ground state condensation energies (b = 0, 1) as a function of d/A for 
A = 0.224.. Q ranges from 22 ( 23) up to 400 (401) for even (odd) grains and m = 60. 
The PBCS results are those of ref. [S7j- Taken from [21]. Reprinted by permission of 
the American Physical Society. 

question of what is the fate of the Josephson effect between two small superconducting 
grains. This question has recently been addressed in ref. [55] . 

To study the physics of two very small grains, it is necessary to have a computational 
framework that can accurately describe the dynamics. From what we discussed earlier, 
BCS approximation is inadequate in a regime in which the two grains are sufficiently 
small. The Richardson approach, which gives the exact solution for individual grains, 
no longer applies when there are two distinct grains with tunnelling between them. 
The phDMRG approach, which had been used successfully for individual grains, can 
be naturally generalized to a two-grain scenario and might therefore provide a useful 
framework for an accurate treatment of the tunnelling phenomenon. In ref. jH], the 
two-grain phDMRG method was developed and applied to this problem. The basic idea 
of the method is to add a particle level and a hole level for one type of grain and then 
to follow that by adding a particle level and a hole level of the other type, gradually 
moving away from the Fermi surface. At each step, truncation is implemented in the 
particle block (which includes levels from both grains) and the hole block (likewise with 
levels from both grains) based on the associated reduced density matrices. This process 
is continued until all levels have been treated. 

There are two variables that define the two-grain problem. One is the size of the 
grains, or equivalently the number of single-particle levels that contribute. The second 
is the Josephson coupling, namely the strength of the term in the Hamiltonian that 
governs pair tunnelling. 

For the purposes of this review, perhaps the key point that emerged from the study 
concerns the range of validity of the two-grain DMRG method as a function of these 
two variables. This is schematically illustrated in fig. [TUl together with an analogous 
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Figure 10. Schematic illustration of the regimes of validity of the two-grain DMRG 
method and the tight-binding (tb) method in the description of tunnelling between 
two superconducting grains. 7 represents the inter-grain coupling and n the number 
of energy levels. Taken from Ref. |69j . 

illustration of the range of validity of the tight-binding (tb) method. 

The tight-binding method - often called the weak-coupling method - assumes that 
the coupled system can be represented as a product of the ground states of two decoupled 
systems, which themselves can be obtained from the ordinary one-grain DMRG method. 
As the figure indicates, this should only work well for very weakly-coupled grains. 

The two-grain DMRG method has a very different range of validity. It does not 
in general work well for weakly-coupled grains, except if the grain sizes are sufficiently 
small. This can be understood from the following argument. The DMRG method, in 
all of its manifestations, depends critically on the presence of correlations between the 
block and the medium in order to effectively truncate the space. When the two grains 
are very weakly coupled, there is no natural way to truncate the block, as all states 
are of comparable importance. Of course, if the number of levels is sufficiently small, 
the number of states m required to adequately describe it is also very small and thus 
the DMRG can still well approximate it. These points are all nicely borne out in the 
calculations reported in ref. [d^] . 

4. 2D Electrons in a High Magnetic Field 

Electrons in high magnetic fields are constrained to move in an effective two-dimensional 
space characterized by the Landau orbits. Depending on the filling of the orbits, various 
phases have been predicted and subsequently observed. The Fractional Quantum Hall 
(FQH) effect, discovered for various fractional fillings in the lowest Landau level, was 
followed by the observation of different exotic phases, like stripes, bubbles, Wigner 
crystallization, etc... 
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An important understanding of these phases has been obtained by exact 
diagonalization for small systems [70] and by Quantum Monte Carlo studies |71j . 
Clearly, there are strong limitations to both methods. Large scale diagonalization can 
only handle systems with N e ~ 10 electrons, while Monte Carlo calculations only provide 
the bulk properties of the ground state, with the accuracy of its results being limited 
by the well-known sign problem. 

The DMRG offers here a great opportunity to go well beyond these limits in the 
study of 2D electrons in the lower Landau levels for different filling factors. As in other 
applications of the DMRG to 2D systems, the main problem is the mapping onto an 
effective ID path. In ref. (2H|, the authors proposed the use of the eigenstates of free 
electrons in a perpendicular magnetic field in the Landau gauge to represent the local 
orbitals that play the role of the single-electron states. The Hamiltonian representing 
the Coulomb repulsion among the electrons in the local orbital basis is written as 

H = S^^c^Cn + — A„ irt2rt3n3 cJ ti cJ l2 c n4 c n3 , (27) 

where S is a classical Coulomb energy and c' n creates an electron in the single-particle 
state n = 1, ...,M, M being the total number of single-particle states in the unit cell, 
and A is the interaction matrix of the Coulomb force (for details see |32j). 

The Hamiltonian (|2*7|) represents a strongly interacting electron system with 
degenerate single-particles states. While the Landau basis provides a natural mapping 
to a ID system by choosing the single-particle states n as the effective sites in a ID 
lattice, the complete degeneracy of the states without the existence of a Fermi energy 
impedes the separation into particle and hole blocks and the identification of a region of 
the single-particle space with strong correlations. Furthermore, the interaction matrix 
A is "long range" in the single-particle space, in the sense that it has finite matrix 
elements between distant sites in the ID lattice. 

The finite algorithm DMRG was applied to the system described by the Coulomb 
Hamiltonian (|2Tj) in the basis of Landau orbits N — 0,1 and 2, to study the ground 
state properties and to determine the phase diagram as a function of the filling factor 

v = N e /M jHBISS]. 

In the warm up, the first left block was constructed from the states n = 1, 2 and 
the first right block from the states n = M, M — 1. The sites with n = 3 and n = M — 2 
were subsequently added to the left and right blocks, respectively, in the first iteration, 
with a truncation based on the eigenvalues of the density matrix implemented to retain 
m states. The method continued in this way until all sites were taken into account, 
completing the warm up. Then, several sweeps were performed until convergence in 
the ground-state energy was obtained. A typical example of the exponential decrease 
of the density matrix eigenvalues is given in fig. El for a system of 18 electrons in 
M = 54 single-particle states. Retaining m = 200 states in each iteration leads to a 
global accuracy of 10~ 4 . 
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Figure 12. The ground-state phase diagram of the third- lowest Landau level obtained 
from a DMRG calculation with up to 25 electrons. Taken from 



To identify the properties of the ground state for the different filling factors, the 
expectation value of the pair correlation function in the DMRG wave function, 

9 fr) = Ne 1} (*l E H* + *i ~ Rj) I*) , (28) 

was used. Here, Rj is the guiding center of the ith electron, r is the relative distance 
between the pair of electrons and L x ^ y are the sides of the unit cell. 

A thorough study of the phase diagram of 2D electrons in the lowest three Landau 
levels was carried out in [SU ESj- We present in fig. H21 as an example, the phase 
diagram of the third Landau level (N = 2) (23] - The diagram shows a Wigner Crystal 
for low filling factors up to v = 0.24, after which a transition to a 2-electron bubble 
phase takes place. There is a second transition to a stripe phase at v — 0.38. The 
properties of the ground state are nicely displayed in fig. [T21 where we can see the quite 
marked geometrical characteristics of the three phases in the pair correlation function. 

The finite algorithm DMRG in the local Landau basis has been shown to constitute 
a powerful tool for investigating the physics of 2D electrons in the lowest Landau levels. 
Further optimization of the method could allow the study of larger systems and the 
investigation of the behavior of the systems close to the phase transitions. 
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Figure 13. Pair correlation function of (a) a stripe state, (b) a two-electron bubble 
state, and (c) a Wigner crystal. Taken from . 



4-1. The DMRG in nuclear structure 

The successes of the DMRG method in Quantum Chemistry, in the treatment of 
ultrasmall superconducting grains and in the study of 2D electronic systems suggests 
that it might also prove useful in the treatment of another finite Fermi system, the 
atomic nucleus. The typical shell-model problem in nuclear structure involves solving 
the Schrodinger equation with a one- and two-body interaction in a spherical shell- 
model basis. The system usually involves two kinds of interacting particles, neutrons 
and protons, with a Hamiltonian 

H = H P + H n + V pn (29) 

that includes three terms - one that acts solely in the proton sector, one that acts solely 
in the neutron sector, and the proton-neutron interaction. 

For heavy enough nuclei and for nuclei with a fairly large number of active (or 
valence) nucleons, exact diagonalization of the shell-model Hamiltonian is not feasible. 
As an example, to treat the J 71 " = 2 + states in the deformed nucleus 154 Sm assuming an 
inert 132 Sn core, 12 valence neutrons in the orbits 3/V/2, 3p 3 / 2 , lh 9 / 2 , U13/2, 3pi/ 2 , 2/ 5 / 2 
and 8 valence neutrons in the orbits 1^7/2, 2rf 5 / 2 , l/in/2, 2rf 3 / 2 , 3si/2 requires 
diagonalization in a basis involving 3.46 x 10 14 states, well beyond the limits of the best 
diagonalization routines. Might it be possible, however, to iteratively obtain accurate 
approximate solutions to the shell-model problem using the DMRG strategy? This 
question has been discussed in several recent works based on the particle-hole algorithm, 
as we now briefly review. 

A characteristic feature of nuclear structure is the presence of several competing 
collective structures, most importantly those associated with the pairing and quadrupole 
fields. In refs. |27^l27j. the phDMRG was applied to a test model involving both pairing 
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Figure 14. DMRG correlation energies as a function of m for a system of 10 particles 
in a j = 25/2 orbit subject to a Hamiltonian with x = lj 9 = 0.1, and e = 0.1. 
The solid line represents an exponential fit to the DMRG results. Taken from |27j . 
Reprinted by permission of the American Physical Society. 

and quadrupole collectivity to assess whether the algorithm is able to handle their 
interplay. For simplicity, the calculations were carried out in a model involving one type 
of nucleon in a large single-j shell, interacting via a P + Q Hamiltonian 

H = -xQ ■ Q ~ gP ] P ~ e H H a )m a im , 

m 

where 

Q» = J2^ m + v\Q»\J m ) c i l+ n c m 

m 

and 

pt = j2(-y- m cid m . (so) 

m 

The third term in the Hamiltonian splits the levels of the single-j shell into a 
set of equally-spaced levels of oblate character, with the largest |m| value lowest. 
Because of that term, the Hamiltonian is not in general rotationally invariant and 
thus its eigenstates do not conserve angular momentum. This unphysical assumption 
was introduced so that each single-particle orbit would be doubly degenerate, greatly 
facilitating the use of an iterative computational algorithm in the test calculation. 

The first results reported were for of 10 particles in a j = 25/2 orbit. 

The calculations were restricted to states with total angular momentum projection 
M — 0, for which the model admits 109,583 states, well within the limits of standard 
diagonalization routines. Thus, for this problem comparison of the DMRG results 
with the exact results was feasible. The results we show were obtained for the case 
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Figure 15. Schematic illustration of the splitting of the spherical single-particle 
levels of the 2s— Id shell into a set of doubly-degenerate levels by an axially-deformed 
Hartree-Fock calculation. The dashed line represents the Fermi energy (Ep), which 
separates the particle levels from the hole levels. Each doubly-degenerate level is 
labelled by its angular momentum projection on the intrinsic z-axis. Taken from ref. 

M 



of x = 1, g — 0.1, e = 0.1, for which Hartree Fock Bogolyubov (HFB) approximation 
shows a well-defined superconducting solution. As such, the test Hamiltonian describes 
a scenario in which pairing and quadrupole correlations indeed compete. 

Figure El compares the correlation energies, defined as the gain in ground-state 
energy relative to Hartree Fock approximation, from the DMRG method and exact 
Lanczos diagonalization. For comparison, the exact and HFB values are E^^ ct = 
-0.70633 and E^ r B = -0.20641. By an m value of 200, the DMRG method produced 
more than 99% of the exact correlation energy, a dramatic improvement over HFB 
approximation which only achieved 28%. The largest superblock matrix that had to be 
diagonalized for this value of m had a dimension of 2, 886. 

Though the calculations used to generate fig. E] only targeted the ground state, 
they produced very accurate results for low-lying excited states as well. For m = 200, 
the lowest three excitation energies were reproduced to 1% or better. 

Reference |2Z1 also reported results for a much larger problem, involving 40 particles 
in a j = 99/2 orbit. In this case, the total number of M = states was 3.84 x 10 25 , 
much too large for exact diagonalization. Here too the results were found to converge 
exponentially, and the estimated error for the ground state correlation energy at m — 100 
was just one part in 10 4 . 

The first reported attempt at a realistic application of the phDMRG method in 
nuclear structure was presented in ref. [2H]- In that work, the nucleus 2A Mg, with four 
neutrons and four protons outside doubly-magic 16 0, was studied. As is usual in the 
nuclear shell model, the Is and Ip core orbits of 16 were assumed to be filled and inert 
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Figure 16. The calculated energies of the ground state and the three lowest excited 
states of 24 Mg as a function of m, the number of states retained in each block. The 
horizontal solid lines refer to the exact energies for the four states. Taken from ref. 
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and the remaining 8 nucleons were restricted to the orbits of the 2s — Id shell. This 
shell-model problem is small and easily solved by exact shell-model diagonalization. 

The results reported in ref. [2Bj were obtained using two different prescriptions for 
the ordering in which single-particle levels were included in the phDMRG algorithm. 
Here we show the results obtained using an axially-symmetric Hartree Fock prescription, 
with the active levels illustrated in fig. ED The axially-symmetric HF single-particle 
levels are all doubly degenerate, thereby facilitating the use of the usual iterative DMRG 
scheme for adding levels. The calculations were done using the so-called USD effective 
Hamiltonian, which has been tailored to nuclei in the 2s — Id shell. Note that the 
DMRG method employed in these calculations is closely related to the two-grain DMRG 
reported in jSH] and discussed in Sect. III.B. 

Figure Uni shows the results for the energies of the lowest four states of the nucleus. 
For this problem the largest m value that can be realized is m = 64. The key point 
to note is that the convergence to the exact result (shown by the solid line) is very 
slow. A value of m above 40 is required to get acceptably accurate reproduction of the 
exact results. Unfortunately, for such values of m, the matrices that must be treated 
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are not much smaller than those of the full problem. Clearly, the phDMRG, method, 
as implemented in ref. [2E], does not work very well for 2A Mg. 

There are several possible reasons for the failure of these calculations. One is that 
the infinite algorithm that was used cannot capture the correlations between the different 
single-particle levels in the problem. Unfortunately, it is difficult to incorporate sweeping 
in the phDMRG, especially for systems with active neutrons and protons, because of 
the preponderance of blocks involved. A second is that the method, as used, does 
not preserve rotational invariance. In the following section, we discuss how rotational 
invariance, and other non-Abelian symmetries, can be built into the DMRG. A third 
concerns the order in which the levels were included. As noted in sect. 13 .1\ a HF 
ordering, or indeed any ordering based solely on single-particle energy considerations, 
may not be optimal. Finally there are possible issues regarding implementation of the 
DMRG truncation algorithm. In particular, no effort was made in the calculation of ref. 

to ensure that at least one state of each partition (number of particles/holes and 
their spin projection) was maintained in each iteration. In a more traditional DMRG 
study of M Mg, it was found that implementation of such a condition dramatically 
improved convergence of the results [72] . 

5. Restoration of Symmetries 

5.1. General remarks 

Symmetries are a crucial element of all quantum systems. The invariance of a system 
under a group of symmetry transformations gives rise to conserved quantum numbers 
for the physical states of the system and associated degeneracies. In the case of uniform 
quantum lattices, translational invariance is a natural symmetry, with the associated 
conserved quantum number being the linear momentum. In spin lattices, the total spin 
of the system is likewise a conserved quantum number for all physical states. And 
in finite Fermi systems, total angular momentum is conserved as a reflection of the 
rotational invariance of the system. 

Methods of approximately solving the Schrodinger equation frequently break the 
symmetries of the Hamiltonian and thereby lead to states that do not contain the 
appropriate quantum numbers. Often this is done intentionally, as spontaneous breaking 
of symmetries is a natural way to include collective correlations between the particles 
in the system. It is then critical, however, to restore the broken symmetries in order to 
establish contact with the states of the physical system, i.e., those probed in experiment. 

The DMRG method, being a truncation of the full Hilbert space of the problem, 
in general violates symmetries, unless it is implemented very carefully. To understand 
this, consider the problem of a spin lattice, where we gradually add spin sites as we grow 
the system. To construct states of good total spin, it is necessary to include all the spin 
projections required by the Clebsch Gordan series. Truncation to a limited number of 
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states based solely on density matrix considerations will not in general keep all the spin 
configurations needed to reconstruct the full spin eigenstates. 

As pointed out by White [Oj, it is possible to preserve spin symmetry, in principle, 
by careful implementation of the density matrix strategy. In particular, as long as 
the density matrix truncation always includes all states from degenerate multiplets, it 
will preserve spin symmetry. However, because of the numerical limitations inherent in 
computational algorithms, it is extremely difficult to guarantee that we do not eliminate 
states that in a perfect calculation were degenerate with those that are kept. Equally 
important, it is not feasible to fix the number of states kept a priori in such a calculation, 
since that number is dictated by the sizes of the degenerate multiplets. 

For the above reasons, several groups have been looking into alternative ways to 
conserve symmetries within the DMRG algorithm. In this section, we briefly review the 
important progress that has been made. 

The first significant progress along these lines was the work of Sierra and Nishino 
[31]. In their work, they showed that Hamiltonians with a continuous symmetry could 
be recast into the form of Interaction Round a Face (IRF) Hamiltonians, familiar in 
statistical mechanics. They then showed how to develop a DMRG algorithm appropriate 
to IRF Hamiltonians. The method was applied to several spin chains with impressive 
success. Examples include the spin-1/2 [31], spin-1 and spin-2 [32] Heisenberg chains 
and the XXZ chain [31]. In all cases considered, the method achieved significantly 
greater accuracy than the usual DMRG algorithm when keeping the same number of 
states. 

Despite the historical importance of the IRF-DMRG algorithm in the proper 
incorporation of symmetries in the DMRG, it is difficult to imagine its use for more 
complex systems with more complicated symmetries. Where major progress was made 
towards the development of a general framework for incorporating symmetries in the 
DMRG was through the work of McCulloch and Gulacsi j3E]- Their work presented a 
formalism that could be applied to any group of symmetry transformations, through 
proper construction of basis states that preserve the symmetry. They refer to this as 
the non-Abelian DMRG method. In the following subsection, we describe in some detail 
how their method can be applied to a finite Fermi system with conserved total angular 
momentum. We will refer to this specific example of their method as the jDMRG 
algorithm. 

5.2. Details on the jDMRG algorithm 

In htis subsection, we develop the jDMRG formalism, showing in some detail how it 
can be implemented for finite Fermi systems so as to produce exact angular momentum 
eigenstates. For simplicity, we restrict the discussion to systems with but one type 
of particle. For systems with two or more types of particles, atomic nuclei being an 
example, some slight modifications are required. 

We begin by stating the problem. We wish to solve the Schrodinger equation for a 
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set of N identical particles subject to a Hamiltonian written in second-quantized form 

as 

. o 



" j{"j"i) u ■ E X hnhu ( a h a l) (VhY 



(31) 



where j = \J1j + 1. 

Note that we are assuming here that the particles under discussion are created 
(annihilated) in the single-particle state {jm} by the operator aj m (a jm ) and that the 
Irreducible Tensor Operator associated with the single-particle annihilation operator is 
defined by dj m = (— y~ m dj_ m . By working in terms of the Irreducible Tensor Operators 
Oj m and a>j m , we are able to couple operators directly and thus properly reflect the scalar 
nature of the Hamiltonian. Note that we use a coupling notation whereby {■■■)m an d 
[■■■]m refer to coupled Irreducible Tensor Operators of rank J and projection M. 

The N particles will be restricted to a finite set of single-particle levels {jm}, as 
depicted schematically in figure ITTk for five orbits. The orbits are shown in ascending 
order according to their single-particle energies ej. 

The jDMRG procedure involves three stages. The first is a preliminary stage, 
in which all necessary information on the individual active single-particle orbits are 
calculated and stored. The second is the "warm-up" stage, familiar to all DMRG 
algorithms. And the third is the likewise familiar "sweeping stage" . 



5.2.1. The preliminary stage. As the procedure will include the effects of all active 
single-particle levels in an angular-momentum coupled representation, we must store 
the following: 

• All single-shell states \ji, n { , c^, Jj >, for the various active orbits {ji}. Here rij 
refers to the number of particles in the state and a refers to any additional quantum 
numbers needed to fully describe the state of total angular momentum Jj. 

• The one-particle Coefficients of Fractional Parentage (CFPs) for all active orbits. 
These are the expansion coefficients required to express an antisymmetrized n- 
particle state as a product of antisymmetrized states of n — 1 particles and a single 
particle. 

• The reduced (or double-bar) matrix elements of all possibly relevant operators for 
those orbits, namely a], [ajaj]^, [aja^, [{a^ja^ and [{aja^} K {djdj} K }° . 

The reduced matrix elements of dj and a] are very simply related to the one-particle 
CFPs. The reduced matrix elements of the more-complicated product operators can be 
obtained from them using well-known formulae [4*5] . 

5.2.2. The warm-up stage. We now discuss the warm-up phase, in which we build a 
first approximation to the dominant structure within groups of active orbits. As we will 
see, most of the formalism we develop for the warm-up phase will also be needed in the 
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Figure 17. Schematic representation of the steps of the the jDMRG method. 

sweep phase. In fig. I17h . we schematically illustrate the "warm-up" phase for the same 
set of five orbits. 

We start off with the first orbit j\, for which all the single-shell matrix elements 
were calculated in the "preliminary" stage. We then add the second orbit, forming an 
enlarged "block". This is schematically illustrated in fig. IT7b . where the enlarged block 
contains orbits ji and j%. 

Next, we must build the reduced matrix elements of all relevant sub-operators of 
the Hamiltonian for the enlarged block. To see how this is done, we consider a generic 
block involving k orbits, \j% ■ ■ ■ jk, n, a, J). The enlarged block basis, in which orbit 
jk+i has been added, is \ji ■ ■ ■ jk, n, a, K; jk+i, Tn, (3, L (J)). The matrix elements of the 
single-particle creation operator a\ in the enlarged block are 

(ji—jk, n',a',K'; m',f3', U (J')\\ a] \\ji-..jk, n, a, K; j k+1 , m, (3, L (J)) 

,K'+L'+j+i 77/ J K J' L \ 



J K i 



5f3j3>5m,m>5LL> ( — ) % J J' 

x 0'i • • -jk, n \ ^'11 a] \\h ■ ■ -Jk, n, a, K) 4 

i X X X I \K'+L+J'+i j j, 
+ Oaa'Onn'OKK' { — ) JJ 



I! J 1 K' 
J L i 



x (j k+l ,m\f3\L'\\a\\\ Jk+1 ,m,f3,L) . (32) 

It involves a reduced matrix element from the previous iteration (which has already 
been calculated) and a reduced matrix element for the orbit that is being added (which 
was calculated in the preliminary stage). Thus, all information is available to calculate 
this reduced matrix element, and indeed all others, in the enlarged block. 

As in all DMRG algorithms, we now wish to choose out of all the states in the 
enlarged block those to keep. Here the strategy is somewhat different depending on 
whether there is only one kind of particle or two. If there are two kinds of particle, 
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one can couple the corresponding states together, build and diagonalize the superblock 
Hamiltonian and then truncate in the sector of interest by using its reduced density 
matrix, assuming each type of particle is the medium for the other. 

When there is only one kind of particle, as in the example under discussion, an 
alternative strategy must be found. A natural one is to use the Wilson RG strategy of 
diagonalizing the Hamiltonian in the enlarged space and keeping the m lowest eigenstates 
for each value of n and J. 

Following truncation, the reduced matrix elements for the truncated enlarged block 
are calculated and stored on hard disk. This is done for all size blocks, j%, ji — > j 2 , 
ji -> J3, h J4, and ji -> j 5 . 

5.2.3. The sweep stage. We now discuss the sweep stage, whereby we sweep back 
and forth through the single-particle levels to improve the description of each block by 
optimally taking into account its correlation with all other levels. 

The sweep stage is schematically illustrated in fig. ITTb for the same five orbits. 
At the point represented by the figure, we are sweeping backwards through the single- 
particle levels and have treated as a block the last two levels, j$ and j'4. We wish to 
enlarge this block to include level j'3. The medium for this enlarged will be levels ji — ► j 2 . 

The sweep stage involves first a calculation of all relevant matrix elements in the 
enlarged block consisting of levels j'3 — > This can be done using exactly the same 
formulae developed for the warm-up stage, with one example given in Eq. ()32j) . 

Once we have calculated all matrix elements in the enlarged block, we construct 
the states of the superblock by coupling those of the enlarged block (orbits j 3 — > j 5 ) to 
those of the medium (orbits j\ — > j 2 ), viz: \ji ■ ■ ■ j 2 ,n,a, K; j 3 ■ ■ ■ j 5 ,m, (3, L (J)). We 
only construct states with the correct total number of particles N = n + m. 

The calculation of the Hamiltonian matrix in the superblock involves a sum of 
terms, each being a product of a reduced matrix element in the enlarged block (just 
calculated) and a reduced matrix element in the medium (output from the previous 
sweep or the warm up). 

As always, the sweep stage is iterated until acceptable convergence has been 
achieved. 

5.3. Applications of the non-Abelian DMRG formalism 

The first application [HE] of the non-Abelian DMRG method was reported by McCulloch 
and Gulacsi for the one-dimensional Hubbard model. This model has the usual 
spin symmetry and also a pseudo-spin symmetry generated by the set of operators 

l + = Ei(-)M T 4 I~ = EiR^iQt , and P = U n n + n n - !)■ Thus the DMRG 
algorithm can be implemented at three levels: (1) by only including the spin projection 
as a conserved quantum number (called U(l)), (2) by including the total spin as a 
conserved quantum number (called SU{2)) and (3) by considering both the total spin 
and the pseudo-spin as good quantum numbers (called SO (A)). There is also an additive 
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Table 5. Comparison of results for three sets of DMRG calculations for the ground 
state of the half- filled 60-site Hubbard model with t = U = 1. The CPU time is in 
seconds per sweep. (1) U(l) x U (1) refers to a calculation in which number conservation 
and conservation of spin projection were imposed; (2) U(l) x 5(7 (2) also includes total 
spin conservation; (3) 5*0(4) in addition includes pseudo-spin conservation. Extracted 
from Table I of ref . . 



Basis 


rn 


D 


E 


(E-E g )/\E g \ 


CPU 


U{1) x [/(l) 


100 


100 


-76.7484986435 


4.2 x 10~ 5 


10 


U(l) x (7(1) 


300 


300 


-76.7516910404 


6.3 x 10~ 7 


110 


U(l) x 517(2) 


100 


226 


-76.7515581914 


2.3 x 10~ 6 


15 


£7(1) x 5(7(2) 


300 


716 


-76.7517389831 


1.1 x 10~ 8 


158 


50(4) 


100 


526 


-76.7517351742 


6.1 x 10~ 8 


18 


50(4) 


300 


1766 


-76.7517398448 


7.9 x 10~ n 


133 



(or Abelian) symmetry of number conservation (likewise (7(1)), which can be included at 
each of the three levels. Results for the three levels of DMRG application are summarized 
in table El for a half- filled 60-site Hubbard lattice with t = U = 1. 

The table includes a quantity D, which in the case of a no n- Abelian symmetry 
is equal to the number of states in the Abelian analysis required to get the same 
results as in the symmetry-adapted analysis. Thus, for example, a calculation in 
which the full SO (4) symmetry is preserved and in which a total of m=100 states are 
maintained in a block would require m = 526 states for comparable accuracy when only 
the z-projection and the particle number are conserved. The table also indicates the 
computational effort required for these various calculations. Keeping the same number 
of states in a calculation with full symmetry restoration as in a calculation with only 
additive symmetries increases the cost in CPU time very little, but yields an enormous 
improvement in accuracy. 

The results shown in Table 03 are for the the S=0 ground state of the system. 
Calculations of higher spin states do not lead to as dramatic an improvement and cost 
savings. 

Subsequently, the method was applied to the two-dimensional t-J Model for lattices 
up to 24x6 [SZ|. Here the purpose was no longer to simply assess the usefulness of 
symmetry restoration in the DMRG approach, but rather to use the power of the method 
to provide useful insight into this model of possible relevance to high-T c superconductors. 
In particular, it attempted to address the important, but controversial, issue of stripe 
formation in the t-J model. 

The analysis made use of the Liang and Pang approach to 2D lattices |49j . which 
unrolls the 2D lattice onto a ID lattice with long-range interactions. The calculations 
included up to m = 1200 states per block. An indication of the accuracy of the method is 
given in Table |H1 where the results obtained for a 16x6 lattice with full spin conservation 
are compared with those having conserved spin projection only [20] ■ As can be readily 
seen, the calculations with full SU(2) spin conservation lead to significantly lower 



The Density Matrix Renormalization Group for finite Fermi systems 



45 



Table 6. Comparison of results of DMRG calculations for the ground state of a 16x6 
t-J system with J = 0.35, t = 1, eight holes and cylindrical boundary conditions. 
The results labelled U(l) only preserve spin projection and are from ref. |5U) . Those 
labelled SU(2) involve full spin conservation. An extrapolated "true" ground state 
energy is also shown. From Table I of ref. [37]. Reprinted by permission of Taylor & 
Francis Ltd., http://www.tandf.co.uk/journals/titles/14786435 



Basis 


m 


E 


U(l) 


1000 


-52.279 


SU(2) 


500 


-52.284 


SU(2) 


800 


-52.463 


SU(2) 


1200 


-52.520 




00 


-52.65 ±0.05 



energies for comparable block sizes. The ground state emerges from this calculation 
in a two-stripe configuration, in agreement with ref. [5T1] . 

More recently, the method was applied to the problem of localized spin ordering 
in Kondo lattice models [73]. By exactly preserving spin symmetry in the problem, the 
authors were able to directly measure the magnetization of the ground state and this 
made it possible for them to confirm the existence of a second ferromagnetic phase in 
the model [?3], for which weak signals had been observed in a previous DMRG study 
that did not fully restore symmetries [7H] . 

6. Summary and Outlook 

In this review, we have described the current status of efforts to develop the Density 
Matrix Renormalization Group for application to finite Fermi system. Enormous 
progress has been made on this front and the method has indeed been applied with 
impressive success in several areas, especially in applications to Quantum Chemistry, to 
study the physics of ultrasmall superconducting grains and to study the properties of 
two-dimensional electron systems. Along with these applications has come much work 
aimed at identifying the key issues that must be resolved in order to build an optimized 
DMRG algorithm for use in these and other systems. One of the most important 
issues concerns the optimal approach to ordering the single-particle basis states that are 
iteratively added by the DMRG algorithm. While important progress has been made, 
questions still remain and more work remains to be done. 

Another area in which important progress has been made concerns the full 
incorporation of symmetries in the DMRG method. It is now possible to consistently 
incorporate non-Abelian symmetries, such as rotational invariance, in addition to 
additive or Abelian symmetries. Early applications of the associated new methodology 
suggest that it permits increased accuracy in DMRG calculations, while requiring less 
CPU cost. More work along these lines is clearly warranted. 
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Another interesting point that has recently emerged concerns the usefulness of 
the DMRG method when treating weakly-coupled strongly-correlated subsystems. 
An application of the particle-hole DMRG method to tunnelling between ultrasmall 
superconducting grains suggests that the methodology breaks down when the tunnel 
coupling between the two grains is very weak. Similar issues also arise when dealing 
with a set of single-particle levels that involve well-separated subsets (e.g. shells), which 
are also very weakly coupled [72]. It would be interesting to see whether the DMRG 
algorithm can be appropriately modified to handle such scenarios. A possibility might 
be to apply the DMRG separately to the various subsystems and then to couple them 
perturbatively. 

The issue of the ordering of levels for their optimal treatment within the DMRG 
method has recently led to the suggestion that Quantum Information Theory might 
provide a useful handle on the problem. This connection needs further exploration. 

We have focussed our remarks on some of the most critical questions that have 
recently been raised in the context of the DMRG method as it applies to finite Fermi 
systems and some of the most interesting methodological breakthroughs. As further 
progress is made to address these and other issues regarding the optimal protocol for 
DMRG studies, we expect a rapid increase in the number of applications. A potentially 
fertile area that has not yet been explored is the atomic nucleus. 

In closing, we believe that the Density Matrix Renormalization Group method will 
continue to grow as a viable method for carrying out reliable and extremely accurate 
calculations for the properties of complex finite Fermi systems, adding to the enormous 
history of success it has already achieved in the description of quantum lattices. 
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